########################################
###           CORRECTION TO:         ###
###           Public Health          ###
###   in the Aftermath of Civil War: ###
###      A Spatial Time-Series       ###
###      Cross-Sectional Analysis    ###
###           July 19, 2025          ###
###           Yuta KAMAHARA          ###
###        Public Choice Studies     ###
########################################

# This R script is prepared for the correction of Figure 1 in my article.
rm(list = ls()); gc()

getwd()

# folders are created.
dir.create("fig")

# required packages 
library("cshapes") # create the contiguity matrix
library("readstata13") # to read Stata13's .dta
library("tidyverse") # `dplyr`; `ggplot2`
library("sf") # st_as_sf for map
library("sp") # sp::plot
library("ggpubr") # multiple maps in one figure.

qog <- read.dta13("https://www.qogdata.pol.gu.se/dataarchive/qog_std_ts_jan19.dta") # Data Archive of The QoG Institute
warnings()

# When using HALE, 
qog.5yr <- qog %>% 
  dplyr::select(cname, year, ccodecow, ccodewb, who_halet) %>%
  dplyr::filter(year>=1951, year<2016)

# Old and new Sudan exist simultaneously. So, I deleted this duplication.
sudan11 <- qog.5yr %>% 
  dplyr::filter(cname=="Sudan (-2011)"&year <= 2011)
sudan12 <- qog.5yr %>% 
  dplyr::filter(cname=="Sudan (2012-)"&year > 2011)
qog.5yrs <- qog.5yr %>% 
  dplyr::filter(cname!="Sudan (-2011)"&cname!="Sudan (2012-)")

qog.5yr <- rbind(qog.5yrs, sudan11, sudan12)

rm(sudan11,sudan12,qog.5yrs)

qog.5yr$cname[qog.5yr$cname=="Sudan (-2011)"] <- "Sudan"
qog.5yr$cname[qog.5yr$cname=="Sudan (2012-)"] <- "Sudan"

# Germany(1991-), West Germany(-1990), East Germany(-1990)
gmy90 <- qog.5yr %>% 
  dplyr::filter(cname=="Germany"&year>1990)

gmye <- qog.5yr %>% 
  dplyr::filter(cname=="Germany, East"&year<1991)

gmyw <- qog.5yr %>% 
  dplyr::filter(cname=="Germany, West"&year<1991)

qog.5yrg <- qog.5yr %>% 
  dplyr::filter(cname!="Germany"&cname!="Germany, West"&cname!="Germany, East")

qog.5yr <- rbind(qog.5yrg, gmye, gmyw, gmy90)
qog.5yr$cname[qog.5yr$cname=="Germany, West"] <- "Germany"

rm(gmy90, gmye, gmyw, qog.5yrg)

# Cyprus (-1974), Cyprus (1975-)
cyp75 <- qog.5yr %>% 
  dplyr::filter(cname=="Cyprus (1975-)"&year>1974)

cyp74 <- qog.5yr %>% 
  dplyr::filter(cname=="Cyprus (-1974)"&year<1975)

qog.5yrc <- qog.5yr %>% 
  dplyr::filter(cname!="Cyprus (1975-)"&cname!="Cyprus (-1974)")

qog.5yr <- rbind(qog.5yrc, cyp74, cyp75)
qog.5yr$cname[qog.5yr$cname=="Cyprus (-1974)"] <- "Cyprus"
qog.5yr$cname[qog.5yr$cname=="Cyprus (1975-)"] <- "Cyprus"

rm(qog.5yrc, cyp74, cyp75)

# Ethiopia (-1992), Ethiopia (1993-)
eth92 <- qog.5yr %>% 
  dplyr::filter(cname=="Ethiopia (-1992)"&year<1993)

eth93 <- qog.5yr %>% 
  dplyr::filter(cname=="Ethiopia (1993-)"&year>1992)

qog.5yre <- qog.5yr %>% 
  dplyr::filter(cname!="Ethiopia (-1992)"&cname!="Ethiopia (1993-)")

qog.5yr <- rbind(qog.5yre, eth92, eth93)
qog.5yr$cname[qog.5yr$cname=="Ethiopia (-1992)"] <- "Ethiopia"
qog.5yr$cname[qog.5yr$cname=="Ethiopia (1993-)"] <- "Ethiopia"

rm(qog.5yre, eth92, eth93)


# France (-1962), France (1963-)
fra62 <- qog.5yr %>% 
  dplyr::filter(cname=="France (-1962)"&year<1963)

fra63 <- qog.5yr %>% 
  dplyr::filter(cname=="France (1963-)"&year>1962)

qog.5yrf <- qog.5yr %>% 
  dplyr::filter(cname!="France (-1962)"&cname!="France (1963-)")

qog.5yr <- rbind(qog.5yrf, fra62, fra63)
qog.5yr$cname[qog.5yr$cname=="France (-1962)"] <- "France"
qog.5yr$cname[qog.5yr$cname=="France (1963-)"] <- "France"

rm(qog.5yrf, fra62, fra63)

# Malaysia (-1965), Malaysia (1966-)
mal65 <- qog.5yr %>% 
  dplyr::filter(cname=="Malaysia (-1965)"&year<1966)

mal66 <- qog.5yr %>% 
  dplyr::filter(cname=="Malaysia (1966-)"&year>1966)

qog.5yrm <- qog.5yr %>% 
  dplyr::filter(cname!="Malaysia (-1965)"&cname!="Malaysia (1966-)")

qog.5yr <- rbind(qog.5yrm, mal65, mal66)
qog.5yr$cname[qog.5yr$cname=="Malaysia (-1965)"] <- "Malaysia"
qog.5yr$cname[qog.5yr$cname=="Malaysia (1966-)"] <- "Malaysia"

rm(qog.5yrm, mal65, mal66)

# Pakistan (-1970), Pakistan (1971-)
pak70 <- qog.5yr %>% 
  dplyr::filter(cname=="Pakistan (-1970)"&year<1971)

pak71 <- qog.5yr %>% 
  dplyr::filter(cname=="Pakistan (1971-)"&year>1970)

qog.5yrp <- qog.5yr %>% 
  dplyr::filter(cname!="Pakistan (-1970)"&cname!="Pakistan (1971-)")

qog.5yr <- rbind(qog.5yrp, pak70, pak71)
qog.5yr$cname[qog.5yr$cname=="Pakistan (-1970)"] <- "Pakistan"
qog.5yr$cname[qog.5yr$cname=="Pakistan (1971-)"] <- "Pakistan"


rm(qog.5yrp, pak70, pak71)

# Vietnam (1977-), Vietnam, North (-1976), Vietnam, South (-1976)
vien <- qog.5yr %>% 
  dplyr::filter(cname=="Vietnam, North"&year<1977)

vies <- qog.5yr %>% 
  dplyr::filter(cname=="Vietnam, South"&year<1977)

viev <- qog.5yr %>% 
  dplyr::filter(cname=="Vietnam"&year>1976)

qog.5yrv <- qog.5yr %>% 
  dplyr::filter(cname!="Vietnam, North"&cname!="Vietnam, South", cname!="Vietnam")

qog.5yr <- rbind(qog.5yrv, vien, vies, viev)
qog.5yr$cname[qog.5yr$cname=="Vietnam, North"] <- "Vietnam"

rm(qog.5yrv, vien, vies, viev)

# Yemen (1990-), North Yemen (-1989), South Yemen
yemn <- qog.5yr %>% 
  dplyr::filter(cname=="Yemen, North"&year<1990)

yems <- qog.5yr %>% 
  dplyr::filter(cname=="Yemen, South"&year<1990)

yemy <- qog.5yr %>% 
  dplyr::filter(cname=="Yemen"&year>1989)

qog.5yry <- qog.5yr %>% 
  dplyr::filter(cname!="Yemen, North"&cname!="Yemen, South", cname!="Yemen")

qog.5yr <- rbind(qog.5yry, yemn, yems, yemy)

rm(qog.5yry, yemn, yems, yemy)


# Serbia == 342, (South Sudan == 626); 
# Tibet == 711; Germany, West == 260; Vietnam, South == 817; Yemen, North == 678; Yemen, South == 680
qog.5yr$ccodecow[qog.5yr$cname=="Sudan"&is.na(qog.5yr$ccodecow)] <- 625
qog.5yr$ccodecow[qog.5yr$cname=="Serbia"] <- 342
qog.5yr$ccodecow[qog.5yr$cname=="Tibet"] <- 711
qog.5yr$ccodecow[qog.5yr$cname=="Cyprus"]<- 352
qog.5yr$ccodecow[qog.5yr$cname=="Ethiopia"] <- 530
qog.5yr$ccodecow[qog.5yr$cname=="Germany"] <- 255 # Original code of Germany, West is 260
qog.5yr$ccodecow[qog.5yr$cname=="France"] <- 220
qog.5yr$ccodecow[qog.5yr$cname=="Malaysia"] <- 820
qog.5yr$ccodecow[qog.5yr$cname=="Pakistan"] <- 770
qog.5yr$ccodecow[qog.5yr$cname=="Vietnam"] <- 816
qog.5yr$ccodecow[qog.5yr$cname=="Vietnam, South"] <- 817
qog.5yr$ccodecow[qog.5yr$cname=="Yemen, North"] <- 678
qog.5yr$ccodecow[qog.5yr$cname=="Yemen, South"] <- 680

# I must recode country name because Cyprus and other countries have duplicated observations

# setup for making new ID
trunc((2018-1951)/5+1) # check
trunc((1951-1951)/5+1)
trunc((1955-1951)/5+1)
trunc((1956-1951)/5+1)
trunc((1959-1951)/5+1)
trunc((1964-1951)/5+1)
trunc((1965-1951)/5+1)
trunc((2015-1951)/5+1)
trunc((2016-1951)/5+1)
qog.5yr <- qog.5yr %>% 
  dplyr::mutate(ID.period=trunc((year-1951)/5+1))
table(qog.5yr$ID.period)

# Y: 1960 (representatives for 1956-1960, t-1 for 1951-1955);
# Y: 1965 (representatives for 1961-1965, t-1 for 1956-1960);
# Y: 1970 (representatives for 1966-1970, t-1 for 1961-1965);
# Y: 1975 (representatives for 1971-1975, t-1 for 1966-1970);
# Y: 1980 (representatives for 1976-1980, t-1 for 1971-1975);
# Y: 1985 (representatives for 1981-1985, t-1 for 1976-1980);
# Y: 1990 (representatives for 1986-1990, t-1 for 1981-1985);
# Y: 1995 (representatives for 1991-1995, t-1 for 1986-1990);
# Y: 2000 (representatives for 1996-2000, t-1 for 1991-1995);
# Y: 2005 (representatives for 2001-2005, t-1 for 1996-2000);
# Y: 2010 (representatives for 2006-2010, t-1 for 2001-2005);
# Y: 2015 (representatives for 2011-2015, t-1 for 2006-2010);

qog.5yr1 <- qog.5yr %>% 
  dplyr::rename(ccodecow2 = ccodecow, ID.period2=ID.period)

qog.5yr1 <- qog.5yr1 %>% 
  aggregate(by=list(ccodecow=qog.5yr1$ccodecow2, ID.period=qog.5yr1$ID.period2), mean, na.rm=TRUE) 

View(qog.5yr1)

summary(qog.5yr1$year)

df <- qog.5yr1 %>% 
  dplyr::select(-cname) %>%  
  dplyr::mutate(lbl.period=ID.period, lbl.period=factor(ID.period, levels = c(1:13),
                                                        labels= c("1951-'55", "1956-'60", "1961-'65", "1966-'70", 
                                                                  "1971-'75","1976-'80", "1981-'85", "1986-'90", "1991-'95",
                                                                  "1996-2000", "2001-'05", "2006-'10", "2011-'15")))

df[is.na(df)] <- NA # change from NaN to NA


###################
#### Figure C.1 ###
###################

### 2000

map2000 <- cshp(date=as.Date("2000-12-31"), useGW = FALSE) # use COWCODE
sp::plot(map2000$geometry)

df2000 <- df %>% 
  dplyr::filter(ID.period==10) %>% 
  dplyr::filter(!is.na(ccodecow))

df2000$ccodecow <- as.numeric(as.character(as.factor(df2000$ccodecow)))

analysis2000 <- map2000 %>% 
  dplyr::left_join(df2000, by = c("cowcode" = "ccodecow"))

exp2000 <- round(mean(analysis2000$who_halet, na.rm = TRUE),4)
border2000 <- cshp(date=as.Date("2000-12-31"))
map2000 <- ggplot(st_as_sf(analysis2000)) + 
  geom_sf(aes(fill=who_halet)) +
  geom_sf(data=border2000, color="black", fill=NA) +
  theme_void() + 
  labs(title=bquote("1996-2000"  ~ "("~ bar(H)[2000] == .(exp2000) ~")")) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_fill_gradient(low='black', high='grey95') +
  guides(fill=guide_legend(title="HALE"))

### 2005

map2005 <- cshp(date=as.Date("2005-12-31"), useGW = FALSE) # use COWCODE
sp::plot(map2005$geometry)

df2005 <- df %>% 
  dplyr::filter(ID.period==11) %>% 
  dplyr::filter(!is.na(ccodecow))

df2005$ccodecow <- as.numeric(as.character(as.factor(df2005$ccodecow)))

analysis2005 <- map2005 %>% 
  dplyr::left_join(df2005, by = c("cowcode" = "ccodecow"))

exp2005 <- round(mean(analysis2005$who_halet, na.rm = TRUE),4)
border2005 <- cshp(date=as.Date("2005-12-31"))
map2005 <- ggplot(st_as_sf(analysis2005)) + 
  geom_sf(aes(fill=who_halet)) +
  geom_sf(data=border2005, color="black", fill=NA) +
  theme_void() + 
  labs(title=bquote("2001-2005"  ~ "("~ bar(H)[2005] == .(exp2005) ~")")) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_fill_gradient(low='black', high='grey95') +
  guides(fill=guide_legend(title="HALE"))

### 2010

map2010 <- cshp(date=as.Date("2010-12-31"), useGW = FALSE) # use COWCODE
sp::plot(map2010$geometry)

df2010 <- df %>% 
  dplyr::filter(ID.period==12) %>% 
  dplyr::filter(!is.na(ccodecow))

df2010$ccodecow <- as.numeric(as.character(as.factor(df2010$ccodecow)))

analysis2010 <- map2010 %>% 
  dplyr::left_join(df2010, by = c("cowcode" = "ccodecow"))

exp2010 <- round(mean(analysis2010$who_halet, na.rm = TRUE),4)
border2010 <- cshp(date=as.Date("2010-12-31"))
map2010 <- ggplot(st_as_sf(analysis2010)) + 
  geom_sf(aes(fill=who_halet)) +
  geom_sf(data=border2010, color="black", fill=NA) +
  theme_void() + 
  labs(title=bquote("2006-2010"  ~ "("~ bar(H)[2010] == .(exp2010) ~")")) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_fill_gradient(low='black', high='grey95') +
  guides(fill=guide_legend(title="HALE"))

### 2015

map2015 <- cshp(date=as.Date("2015-12-31"), useGW = FALSE) # use COWCODE
sp::plot(map2015$geometry)

df2015 <- df %>% 
  dplyr::filter(ID.period==13) %>% 
  dplyr::filter(!is.na(ccodecow))

df2015$ccodecow <- as.numeric(as.character(as.factor(df2015$ccodecow)))

analysis2015 <- map2015 %>% 
  dplyr::left_join(df2015, by = c("cowcode" = "ccodecow"))

exp2015 <- round(mean(analysis2015$who_halet, na.rm = TRUE),4)
border2015 <- cshp(date=as.Date("2015-12-31"))
map2015 <- ggplot(st_as_sf(analysis2015)) + 
  geom_sf(aes(fill=who_halet)) +
  geom_sf(data=border2015, color="black", fill=NA) +
  theme_void() + 
  labs(title=bquote("2011-2015"  ~ "("~ bar(H)[2015] == .(exp2015) ~")")) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_fill_gradient(low='black', high='grey95') +
  guides(fill=guide_legend(title="HALE"))

ggpubr::ggarrange(map2000, map2005, map2010, map2015, ncol=2, nrow=2,  
                  common.legend = TRUE,legend="bottom")

ggsave("fig/figC1.png", width=153, height=116,  units = 'mm', bg = 'white', dpi = 300)

dev.off()
